#include "manager.h"
#include "simulatorParameters.h"

void computeDij(pFace face, pGEntity gEntity);

// Calculates Cij, Dij(boundary faces) and nodal volume
// It also works for multi domains meshes
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
void calculateCoefficients3D_old(pMesh theMesh, set<int> &setOfDomains, const char *meshFilename)
{
	cout << "calculateCoefficients3D_old\n";
	cout << "Mesh dimension: 3D\n";
	
	SimulatorParameters simPar;
	simPar.loadParameters();

	vector<pVertex> tetraVertex, edgeVertex;
	int i,j,k,K;

	double tCenter[3], eCenter[3], Cij[3], v[3], *vec[2], val, *tetraFCenter[4];
	double normal[3], dp, I[3], J[3], IJ[3], proj[3];

	pGEntity gEntity;
	pEntity tetra, edge;
	// allocate vectors
	for (i=0; i<4; i++) tetraFCenter[i] = new double[3];
	for (i=0; i<2; i++) vec[i] = new double[3];	

	char *Cij_string[3];
	for (i=0; i<3; i++) Cij_string[i] = new char[64];

	// loop over elements
	// for each element:	1 - calculate Cij for each edge and store them on the respective edge
	// 						2 - calculate element contribution to volume of control volume	

	RIter rit = M_regionIter(theMesh);
	while (tetra = RIter_next(rit))
	{
		// get all four tetrahedron's vertices
		M_GetVertices(tetra,tetraVertex);
		int tetraFaces[4][3] = { {EN_id(tetraVertex[0]),EN_id(tetraVertex[1]),EN_id(tetraVertex[2])},
				{EN_id(tetraVertex[0]),EN_id(tetraVertex[1]),EN_id(tetraVertex[3])},
				{EN_id(tetraVertex[0]),EN_id(tetraVertex[2]),EN_id(tetraVertex[3])},
				{EN_id(tetraVertex[1]),EN_id(tetraVertex[2]),EN_id(tetraVertex[3])}};

		getFCenter(tetraVertex[0], tetraVertex[1], tetraVertex[2], tetraFCenter[0]);
		getFCenter(tetraVertex[0], tetraVertex[1], tetraVertex[3], tetraFCenter[1]);
		getFCenter(tetraVertex[0], tetraVertex[2], tetraVertex[3], tetraFCenter[2]);
		getFCenter(tetraVertex[1], tetraVertex[2], tetraVertex[3], tetraFCenter[3]);

		// get tetrahedron center
		tCenter[0] = tCenter[1] = tCenter[2] = .0;
		R_center(tetra, tCenter);

		// step #1: Cij
		// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		// identify which domain the current tetrahedron belongs to
		// edges on domain's partition have differents Cij, one for each domain 
		gEntity =(pGEntity)R_whatIn(tetra);
		
		const int dom = GEN_tag(gEntity);

		setOfDomains.insert(dom); // store every new domain

		sprintf(Cij_string[0],"Cx-%d",dom);
		sprintf(Cij_string[1],"Cy-%d",dom);
		sprintf(Cij_string[2],"Cz-%d",dom);

		void *r_iter = 0;

		// takes a list of edges for tetra
		pPList edgeList =  R_edges(tetra,0);

		// loop over edges belonging to tetra
		while (edge = (pEntity)PList_next(edgeList,&r_iter) )
		{
			M_GetVertices(edge,edgeVertex);

			Cij[0] = Cij[1] = Cij[2] = .0;			

			// edge's Cij is the sum of two vectors. Both are orthogonals to planes
			// defined by three vectors. They are: 
			// v - ec->tc:		edge center to tetra center
			// vec[0] - ec->fcr:	edge center to right face center
			// vec[1] - ec->fcl:	edge center to left face center
			// Cij = (v)x(vec[0]) + (v)x(vec[1]), Cij points to outside of control volume

			// create vector v:
			for (i=0; i<3; i++) I[i] = J[i]= .0;
			V_coord(edgeVertex[0],I);
			V_coord(edgeVertex[1],J);
			makeVector(I,J,IJ);
			if ( EN_id(edgeVertex[0]) > EN_id(edgeVertex[1]) ) makeVector(J,I,IJ);

			// get edge center
			eCenter[0] = eCenter[1] = eCenter[2] = .0;
			for (i=0; i<3; i++) eCenter[i] = (I[i]+J[i])/2.0;
			//E_center(edge,eCenter);

			// create vector v:
			for (i=0; i<3; i++) v[i] = tCenter[i] - eCenter[i];			

			// search for tetra faces that share the same edge and get their centers.
			// Of course, there are only two faces sharing the same edge.
			// FMDB can do this job easily, but for 3D big meshes a face structure
			// has a high memory cost.

			//int keep;
			K = 0;
			for (i=0; i<4; i++)				// loop over tetra's faces
			{
				for (j=0; j<3; j++)			// loop over face's vertices
				{
					if (EN_id(edgeVertex[0]) == tetraFaces[i][j])
					{
						for (k=0; k<3; k++)
							if (EN_id(edgeVertex[1])==tetraFaces[i][k])
								makeVector(eCenter, tetraFCenter[i], vec[K++]);
					}
				}
			}

			// IJ vector is a reference to Cij. IJ points from node I to node J
			// where: I_id < J_id
			// vec projections on edge IJ must have the same orientation as IJ vector

			double n = norm(IJ);

			// Cij calculation
			for (i=0; i<2; i++)
			{
				for (j=0; j<3; j++) normal[j] = proj[j] = .0;
				cross(v,vec[i],normal);
				val = dot(normal,IJ)/(n*n);
				for (j=0; j<3; j++) proj[j] = val*IJ[j];
				double sinal = ( dot(proj,IJ) < .0)?-1.:1.;
				for (j=0; j<3; j++) Cij[j] += sinal*normal[j]/2.0;
			}			
			// update edge Cij
			for (i=0; i<3; i++)
			{
				val = .0;
				EN_getDataDbl(edge,MD_lookupMeshDataId(Cij_string[i]),&val);
				val += Cij[i];
				EN_attachDataDbl(edge,MD_lookupMeshDataId(Cij_string[i]),val);
			}
			edgeVertex.clear();
		}
		PList_delete(edgeList);

		// step #2: volume of control volume
		// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -		
		const double porosity = simPar.getPorosity(dom);
		
		const double volume = .25*R_Volume(tetra);		 
		char volumeString[64];
		sprintf(volumeString,"volume-%d",dom);
		
		for (i=0; i<4; i++)
		{
			// total nodal volume
			val = .0;
			EN_getDataDbl(tetraVertex[i],MD_lookupMeshDataId(volumeString),&val);
			val += volume;
			EN_attachDataDbl(tetraVertex[i],MD_lookupMeshDataId(volumeString),val);
			
			// weighted porosity
			val = .0;
			EN_getDataDbl(tetraVertex[i],MD_lookupMeshDataId("wght-porosity"),&val);
			val += volume*porosity;
			EN_attachDataDbl(tetraVertex[i],MD_lookupMeshDataId("wght-porosity"),val);
		}
		tetraVertex.clear();
	}
	RIter_delete(rit);	// END TETRAHEDRALS LOOP

	// deallocate vectors
	for (i=0; i<4; i++) delete[] tetraFCenter[i];
	for (i=0; i<2; i++) delete[] vec[i];

	// Dij coefficients are associated to boundary faces.
	// The following loop is made on boundary faces provided by mesh generator.
	// It's expected that all boundary faces be present on mesh file and not 
	// only those with flags. If such faces were not provided, FMDB will create   
	// a data structure to generate them. Note that internal and external faces  
	// to domain will be created. This operation is memory and cpu costly.

	printf("Num. faces: %d\n",M_numFaces(theMesh));
	pEntity face;	
	FIter fit;
	if ( M_numFaces(theMesh) != 0 )
	{		
		fit = M_faceIter(theMesh);
		while (face = FIter_next(fit))
		{
			gEntity =(pGEntity)F_whatIn(face);
			computeDij(face,gEntity);		
		}
	}
	FIter_delete(fit);

	for (i=0; i<3; i++) delete[] Cij_string[i];
	
#ifdef _HAVE_PETSC_
	unifyCijAmongProcessors(theMesh,setOfDomains);
	unifyVolumesAmongProcessors(theMesh,setOfDomains,"volume");
	unifyVolumesAmongProcessors(theMesh,setOfDomains,"wght-porosity");
#endif

	
	string str(meshFilename);	
	string::size_type start = str.find_last_of("/");
	start = (start == string::npos)?0:start;	
	string::size_type end = str.find_last_of(".");

	char buffer[256];
	memset( buffer, '\0', 256 );	
	str.copy(buffer,end-start,start);
	cout<<"buffer: " << buffer << endl;

	char filename[256];
	sprintf(filename,"preprocessor-datafiles/%s-%d_of_%d.dat",buffer,P_pid(),P_size());
	cout<<filename<<endl;
	exportCoefficients(theMesh,3,setOfDomains,filename);

#ifdef _PREPROCESSOR_DEBUG_	
	verifyGeometricCoefficients(theMesh,setOfDomains);
#endif

}
double sum[3]={.0,.0,.0};

void computeDij(pFace face, pGEntity gEntity)
{

	int i;
	char *Dij_string[3];
	for (i=0; i<3; i++) Dij_string[i] = new char[64];	

	sprintf(Dij_string[0],"Dx-%d",GEN_tag(gEntity));
	sprintf(Dij_string[1],"Dy-%d",GEN_tag(gEntity));
	sprintf(Dij_string[2],"Dz-%d",GEN_tag(gEntity));

	pVertex oppositeVertex;	
	findVertexOpposite(face, oppositeVertex);

	double Dij[3] = {.0, .0, .0};
	DijVector(face, oppositeVertex, Dij);

	for (i=0; i<3; i++)
	{
		//printf("Dij[] = {%f %f %f}\n",Dij[0],Dij[1],Dij[2]);
		EN_attachDataDbl(face,MD_lookupMeshDataId(Dij_string[i]),Dij[i]);
		delete[] Dij_string[i];
	}
}

#ifdef _PREPROCESSOR_DEBUG_

// sum(Cij) + sum(Dij) = 0 per domain
void verifyGeometricCoefficients(pMesh theMesh, set<int> &setOfDomains)
{
	int i,j;

	char *Cij_string[3], *Dij_string[3];
	for (i=0; i<3; i++) 
	{
		Cij_string[i] = new char[64];
		Dij_string[i] = new char[64];
	}
	char *Cij_str[3] = {"Cx","Cy","Cz"};
	char *Dij_str[3] = {"Dx-","Dy-","Dz-"};

	int n = 125;

	pEntity edge, face;

	// =========================================================================
	// loop over domains
	set<int>::iterator SIter = setOfDomains.begin();
	for (; SIter != setOfDomains.end(); SIter++)
	{
		double sumcij[3]={.0,.0,.0}, sumdij[3]={.0,.0,.0};
		int dom = *SIter;					// domain flag
		printf("domain %d:\n",dom);

		char str[256];
		double sum=0;
		VIter vit = M_vertexIter(theMesh);
		while (pEntity v = VIter_next(vit))
		{
			double val = .0;
			sprintf(str,"volume-%d",dom);					
			EN_getDataDbl(v,MD_lookupMeshDataId(str),&val);
			sum += val;
		}
		VIter_delete(vit);
		printf("volume dom %d: %f\n",dom,sum);

		// loop over edges domains
		// =====================================================================
		EIter eit = M_edgeIter(theMesh);
		while (edge = EIter_next(eit))
		{			
			vector<pVertex> vertex;
			M_GetVertices(edge,vertex);
			int id0 = EN_id(vertex[0]), id1 = EN_id(vertex[1]);

			if (id0==n || id1==n)
			{
				int id = (id0<id1)?id0:id1;
				double signal = (n > id)?-1.0:1.0;
				// update edge Cij
				for (i=0; i<3; i++)
				{
					double val = .0;
					sprintf(Cij_string[i],"%s-%d",Cij_str[i],dom);					
					EN_getDataDbl(edge,MD_lookupMeshDataId(Cij_string[i]),&val);
					printf("%f ",signal*val);
					sumcij[i] += signal*val;
				}
				printf("\tsum(Cij) = %f %f %f\n",sumcij[0],sumcij[1],sumcij[2]);
			}
			vertex.clear();

		}// end edge loop
		EIter_delete(eit);

		// loop over faces domains
		// =====================================================================
		pGEntity gEntity;
		vector<pEntity> tetras;
		FIter fit = M_faceIter(theMesh);
		while (face = FIter_next(fit))
		{		
			gEntity = (pGEntity)F_whatIn(face);
			int faceFlag = GEN_tag(gEntity);

			getTetraWithFace(tetras,face);
			// if face is external to domain
			if ( faceFlag == 1000 )
			{				
				gEntity = (pGEntity)R_whatIn(tetras[0]);
				int tetraFlag = GEN_tag(gEntity);
				//printf("tetra-flag: %d\n",tetraFlag);
				// check if element belongs to domain
				if (tetraFlag == dom)
				{
					for (i=0; i<3; i++)
					{
						double val = .0;
						sprintf(Dij_string[i],"%s%d",Dij_str[i],faceFlag);
						EN_getDataDbl(face,MD_lookupMeshDataId(Dij_string[i]),&val);
						//sumdij[i] += val;
						printf("%f ",val);
					}	printf("\n");
				}				
			}
			else
			{
				vector<pVertex> vertex;
				M_GetVertices(face,vertex);
				if (EN_id(vertex[0])==n || EN_id(vertex[1])==n || EN_id(vertex[2])==n)
				{
					printf("%d %d %d\n",EN_id(vertex[0]),EN_id(vertex[1]),EN_id(vertex[2]));
					for (j=0; j<tetras.size(); j++)
					{
						gEntity = (pGEntity)R_whatIn(tetras[j]);
						int tetraFlag = GEN_tag(gEntity);
						//	printf("tetra-flag: %d\n",tetraFlag);
						// check if element belongs to domain
						if (tetraFlag == dom)
						{
							pEntity node;
							getVertexOppositeToFace(tetras[j],face,node);
							//printf("\noposite node: %d\n",EN_id(node));
							double tmp[3]={.0,.0,.0};
							for (i=0; i<3; i++)
							{
								sprintf(Dij_string[i],"%s%d",Dij_str[i],faceFlag);
								EN_getDataDbl(face,MD_lookupMeshDataId(Dij_string[i]),&tmp[i]);							
								//printf("%f ",tmp[i]);
							}	//printf("\n");
							double center[3]={.0,.0,.0};
							getFCenter(face,center);
							double nodeCoords[3]={.0,.0,.0};
							V_coord(node,nodeCoords);
							double vec[3] = {center[0]-nodeCoords[0],center[1]-nodeCoords[1],center[2]-nodeCoords[2]};
							double inner_product = .0;
							for (i=0; i<3; i++) inner_product += vec[i]*tmp[i];
							if (inner_product<0) 
								for (i=0; i<3; i++) tmp[i] = -tmp[i];

							for (i=0; i<3; i++) sumdij[i] += tmp[i];
							printf("tmp[]: %.5f  %.5f %.5f\t",tmp[0],tmp[1],tmp[2]);
							printf("sum(Dij) = %.5f  %.5f %.5f\n",sumdij[0],sumdij[1],sumdij[2]);
						}
					}
				}
			}
			tetras.clear();
		}// end face loop
		FIter_delete(fit);
	}	
}

#endif
